Binary Neutron Stars in General Relativity: Quasi-Equilibrium Models 
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We perform fully relativistic calculations of binary neutron stars in quasi-equilibrium circular orbits. 
We integrate Einstein's equations together with the relativistic equation of hydrostatic equilibrium 
to solve the initial value problem for equal-mass binaries of arbitrary separation. We construct 
sequences of constant rest mass and identify the innermost stable circular orbit and its angular 
velocity. We find that the quasi-equilibrium maximum allowed mass of a neutron star in a close 
binary is slightly larger than in isolation. 
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The two-body problem is one of the outstanding, un- 
solved problems in classical general relativity. And yet, 
neutron star binary systems are known to exist, even 
within our own galaxy El. For some of these systems 
(including PSR B1913+16, B1534+12 and B2127-11C) 
general relativistic orbital effects have been measured to 
high precision S. Binary neutron stars are among the 
most promising sources for gravitational wave detectors 
now under construction, like LIGO, VIRGO and GEO. 
This has triggered an intense theoretical effort to predict 
the gravitational wave form emitted during the inspiral 
and coalescence of the two stars. 

Much of the work on binary neutron stars has been per- 
formed within the framework of Newtonian hydrodynam- 
ics H . Many investigators have also studied the problem 
in post-Newtonian (PN) theory. As long as the PN stars 
are well separated, they can be approximated by point 
sources [B, but for close binaries, hydrodynamical effects 
must also be taken into account 

Fully general relativistic treatments of the problem 
are complicated by the nonlinearity of Einstein's equa- 
tions and the requirement of very large computational 
resources. Numerical simulations are currently only in 
their infancy M. Recently, Wilson and Mathews fir] ] 
reported preliminary results obtained with a relativistic 
numerical evolution code. Their dynamical calculations 
suggest that the neutron stars may collapse to black holes 
prior to merger. They also find that, typically, binaries 
have too large a total angular momentum to form a Kerr 
black hole immediately upon merger (see also jll]]). Their 
results are in disagreement with predictions of Newto- 
nian |l^] and PN calculations 0, which show that tidal 
fields stabilize neutron stars against radial collapse. 

In this Letter we report the first calculations in full 
relativity of quasi-equilibrium, equal mass, neutron star 
binaries in synchronized circular orbits. We numerically 
integrate a subset of the Einstein equations, coupled to 
the equations of relativistic hydrodynamics, to solve the 



initial value problem for binaries. Such quasi-equilibrium 
models provide initial data for future dynamical evo- 
lution calculations. We construct quasi-equilibrium se- 
quences of constant rest mass configurations at vary- 
ing separation. These sequences mimic evolutionary se- 
quences in which the stars undergo slow inspiral on nearly 
circular orbits due to the emission of gravitational waves. 
We identify the innermost stable circular orbit (ISCO), 
its angular velocity, and the maximum quasi-equilibrium 
mass of a neutron star in a close binary. 

In Newtonian gravity, strict equilibrium for two stars 
in synchronized circular orbit exists. Since this solu- 
tion is stationary, the hydrodynamical equations reduce 
to the Bernoulli equation, which greatly simplifies the 
problem. Because of the emission of gravitational waves, 
a binary in general relativity cannot be in strict equilib- 
rium. However, outside the ISCO, the timescale for or- 
bital decay by radiation is much longer than the orbital 
period, so that the binary can be considered to be in 
"quasi-equilibrium" . This fact allows us to neglect both 
gravitational waves and wave-induced deviations from a 
circular orbit to very good approximation |[3]]. Some of 
our approximations have been used and calibrated else- 
where Jl4|,[l5| , and a more detailed discussion will be pre- 
sented in a forthcoming paper |fT6f . Here we will briefly 
outline our method and present some of our key results. 

We attempt to minimize the gravitational wave content 
in the solution, in compliance with physical expectations, 
by choosing the 3- metric to be conformally flat [[[o]jl4|] . In 
cartesian coordinates the line element can then be writ- 
ten 

ds 2 = -a 2 dt 2 + *%((fa 4 - uj l dt)(dx J - uj j dt), (1) 

where a is the lapse, cu l the shift and \I/ the conformal 
factor. We satisfy the initial value equations of relativity 
precisely. Our approximation lies in assuming that the 
metric will remain conformally flat for all times. The 
extrinsic curvature Kij then has to satisfy 
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where K 1 ^ — ^ la K l: > and where we have also used the 
maximal slicing condition K l , t = 0. Here V 1 is the flat 
space derivative operator in cartesian coordinates. 

We assume that the matter obeys a polytropic equa- 
tion of state 



P = Kpl +1/n , 



(3) 



where P is the pressure, po the rest-mass density, K the 
polytropic constant, and n the polytropic index. We 
assume that we can neglect deviations from a strictly 
periodic circular orbit and that the stars are corotat- 
ing, which is equivalent to assuming that the fluid four- 
velocity is proportional to a Killing vector. In this 
case the matter equations can be integrated analytically, 
which yields the relativistic Bernoulli equation 
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where q = P/ po, C is a constant of integration and v is 
the proper velocity of the matter. 

The Hamiltonian constraint can now be written 
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Requiring that the maximal slicing condition be main- 
tained at all times, we can use the time evolution equa- 
tion for Kij to find an equation for the lapse, 
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where a = ^a. The momentum constraint becomes 
1 
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where is the constant angular velocity and £ l is a three- 
vector tangent to the matter velocity. With the stars cen- 
tered along the z-axis and orbiting around the j/-axis, we 
have C = (z,0, —x). The last equation can be simplified 
by writing of = G l — jV l B. 

Our approximations reduce the Einstein field equa- 
tions to a set of coupled, quasi-linear elliptic equations 
for the lapse, shift, and the conformal factor (eqs. 
which have to be solved together with the matter equa- 
tion (0) . For boundary conditions at large radius we im- 
pose asymptotic flatness. Solving these equations yields 
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FIG. 1. Rest-density contours in the equatorial plane for a 
neutron star binary close to the ISCO. Each star has a rest 
mass of Mo = 0.178, only slightly below the maximum mass 
at infinite separation, M™ ax = 0.180. The contours span 
densities between the central density and 1% of that value by 
decreasing factors of 0.63. 



a valid solution to the initial value (constraint) equa- 
tions. Such a solution will also provide an approximate 
instantaneous snapshot of a binary evolved according to 
the full Einstein equations, prior to plunge. In the New- 
tonian limit, the above equations reduce to the coupled 
Poisson and Bernoulli equations. 

Our numerical implementation will be described in de- 
tail in Jl(|. Since the stars have equal mass, it is sufficient 
to work in one octant only. We use a full approxima- 
tion storage multigrid scheme to solve the elliptic field 
equations (^-Q) for a given matter distribution. Once 
a solution has been found, the matter can be updated 
(eq. (|J)). This iteration can be repeated until conver- 
gence is achieved to a desired accuracy. We have im- 
plemented this algorithm in a parallel environment using 
DAGH software |L7|]. Typical runs were computed on a 
grid of 64 3 gridpoints. We adjusted the outer boundaries 
for each separation so that the matter was always covered 
by 17 gridpoints along the diameter. 

We determine the rest (baryon) mass Mo, the total 
mass-energy (ADM mass) M, as well as the angular mo- 
mentum J, which refer to the parameters of one individ- 
ual star. Note that physical dimensions enter the prob- 
lem only through the polytropic constant K in (|3|). It is 
therefore convenient to introduce the dimensionlcss quan- 
tities pa = K n p„, M = K- n / 2 M and M = R- n / 2 M. 

In the following we will discuss results for n = 1. A 
survey of several different polytropic indices between 1 
and 2.9 will be presented in |p^| . 

In Fig. 1 we show density profiles for highly relativistic 
neutron stars of rest mass Mo = 0.178 close to the ISCO. 
The maximum mass for such a star in isolation is M ( ? lax = 
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FIG. 2. Rest mass Mo versus central density p c for separa- 
tions za = 0.3 (bottom solid line), 0.2, 0.1 and 0.0 (top line). 
The dashed line is the Oppenheimer-Volkoff result. The insert 
is a blow-up of the region around the maximum masses. 

0.180. Note that the stars at the ISCO are only very 
slightly distorted. In the following we will define the 
ratio between the inner and outer coordinate separation, 

%A ^ in I ^ out- 

In Fig. 2 we plot the allowed rest mass versus the 
central density for several different separations between 
Za = 0.3 (roughly two stellar radii apart) to za = 
(touching). As za — * 1 we expect these curves to ap- 
proach the spherical Oppenheimer-Volkoff (OV) result, 
which we included as the dashed line in Fig. 2. All our 
graphs lie within 2% of the OV curve, showing that the 
presence of a companion star has only very little influence 
on the mass-density equilibrium relationship. 

As we decrease the separation, the mass supported by 
a given p c increases slightly. In particular, the maxi- 
mum mass increases from M™ ax = 0.179 for za = 0.3 to 
Afmax = 0.182 for touching stars. This trend clearly sug- 
gests that the maximum allowed mass of neutron stars 
in close binaries is slightly larger than in isolation. This 
increase is caused partly by the rotation of the stars and 
partly by the tidal fields [fT8[ , Note, however, that we are 
only constructing quasi- equilibrium configurations, which 
may or may not be dynamically stable. For nonrotat- 
ing, isolated stars the maximum mass configuration in 
Fig. 2 marks the onset of radial instability. No gen- 
eral theorem can be trivially applied to binary stars. If 
the results ]l(J are correct (but see |^,|l2|]), it is because 
the equilibrium configurations are dynamically unstable, 
and not because the maximum allowed mass decreases. 
Fig. 2 also shows that keeping the rest mass fixed, the 
central density slightly decreases as the stars approach 
each other and become tidally deformed. Both effects 
are consistent with simple PN predictions Jt]||. 

We construct sequences of constant rest mass M , 
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FIG. 3. The binding energy as a function of the angular 
velocity for several different values of the rest mass Mo. The 
curves are labeled by the compaction (M/R)oo of the stars in 
isolation at infinity. The maximum compaction for a stable, 
isolated, nonrotating n = 1 polytrope is 0.217. The upper 
label gives the frequency for stars with a rest mass of 1.5Mq. 
The dashed curves are the corresponding results from a New- 
tonian version of our code. In the insert we plot results for a 
nearly Newtonian configuration with (M/K)oo = 0.025. Here 
we have also included the results from a Newtonian ellipsoidal 
treatment [12] (long dashed line) and for two point particles 
(dotted line). 

which approximate evolutionary sequences up to the 
ISCO pj| . Note that we have assumed the stars to be 
corotating. This may not be realistic, since it would re- 
quire excessive viscosity pl| . It is more likely that circu- 
lation of the stars is conserved during inspiral, and the 
stars remain nearly irrotational. Nevertheless, we expect 
that our sequences are a reasonable approximation to the 
inspiral up to the ISCO and correctly reveal the effects 
of nonlinear gravitation. 

In Fig. 3 we show the binary binding energy versus an- 
gular velocity for several different rest masses Mq. As the 
stars approach, both finite size effects and nonlinear grav- 
itation play an increasingly important role and cause, for 
stiff enough equations of state, the binding energy to go 
through a minimum and increase again. The location 
of the minimum marks the onset of a secular instability, 
beyond which the binary can no longer maintain circular 
equilibrium. It is expected that the dynamical instability, 
which defines the true ISCO for plunge, occurs beyond, 
but close to, the onset of the secular instability ]l2|] . 

In the Newtonian regime our results agree very well 
with both a Newtonian version of our code and results 
from an ellipsoidal treatment of the binaries [^2| (see in- 
sert). For more relativistic models, comparisons are made 
somewhat ambiguous by the adopted choice of a parame- 
ter to characterize the sequence. Identifying the member 



3 



M 


Moo 


(M/i?)oo 


MqQi sco 


( J tot i J <?(70 


0.112 


0.106 


0.1 


0.01 


1.22 


0.134 


0.126 


0.125 


0.015 


1.12 


0.153 


0.142 


0.15 


0.02 


1.05 


0.169 


0.155 


0.175 


0.025 


1.00 


0.1781 


0.1623 


0.2 


0.03 


0.97 



TABLE I. Numerical values for sequences of constant rest 
mass Mo and polytropic index n = 1. We tabulate the total 
energy M^o and compaction (M/R) oo each star would have in 
isolation at infinity as well as the angular velocity Mof2 and 
the angular momentum J/M 2 at the ISCO. 



at infinity by its value of M (or Mo) versus (M/R)^ (or 
(Mo /R)oo) leads to different Newtonian models and bind- 
ing energy curves, and the differences increase as the stars 
become more compact. In Fig. 3 we choose (M /i?)oo 
and find that the ISCO frequencies agree closely with 
the Newtonian values, but the binding energies differ as 
the compaction increases [Q . 

We summarize these results in Table 1, where we 
also tabulate the dimensionless total angular momen- 
tum J tot /M? ot = J/2M 2 at the ISCO. Note that for high 
enough rest masses this value drops below unity, so that 
these two stars could plunge and form a Kerr black hole 
without having to radiate additional angular momentum. 
Because the orbit will decay rapidly inside the ISCO, its 
presence will leave a measurable imprint on the emitted 
gravitational wave form. Measuring flisco may be the 
crucial ingredient in determining the radius of the star, 
assuming that the mass has been determined during the 
prior inspiral phase |2^] . 
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